home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Mac Magazin/MacEasy 21
/
Mac Magazin and MacEasy Magazine CD - Issue 21.iso
/
Wissenschaft & Technik
/
yorick12vr1-nofpu folder
/
include
/
romberg.i
< prev
next >
Wrap
Text File
|
1995-07-26
|
4KB
|
121 lines
/*
ROMBERG.I
Romberg integrator, after qromb in Numerical Recipes (Press, et.al.)
$Id$
*/
/* Copyright (c) 1994. The Regents of the University of California.
All rights reserved. */
func romberg(function, a, b, epsilon, notvector=)
/* DOCUMENT integral= romberg(function, a, b)
or integral= romberg(function, a, b, epsilon)
returns the integral of FUNCTION(x) from A to B. If EPSILON is
given, Simpson's rule is refined until that fractional accuracy
is obtained. EPSILON defaults to 1.e-6.
If the notvector= keyword is supplied and non-zero, then FUNCTION
may not be called with a list of x values to return a list of
results. By default, FUNCTION is assumed to be a vector function.
If the function is not very smooth, simpson may work better.
SEE ALSO: simpson, max_doublings
*/
{
if (!epsilon || epsilon<0.) epsilon= 1.e-6;
a= double(a);
b= double(b);
s= array(0.0, 5);
h= [1.0, 0.25, 0.0625, 0.015625, 0.00390625];
for (i=1 ; i<=max_doublings ; ++i) {
ss= trapezoid(function, a, b, i, notvector);
if (i >= 5) {
s(5)= ss;
c= d= s; /* Neville algorithm tableau */
ns= 4; /* one less than smallest h, always last */
for (m=1 ; m<5 ; ++m) {
m5= 5-m;
ho= h(1:m5);
hp= h(m+1:5);
w= (c(2:m5+1)-d(1:m5))/(ho-hp);
d(1:m5)= hp*w;
c(1:m5)= ho*w;
dss= (2*ns<m5)? c(ns+1) : d(ns--);
ss+= dss;
}
if (abs(dss) < epsilon*abs(ss)) return ss;
/* extrapolation to h=0 always uses last 5 points */
s(1:4)= s(2:5);
h*= 0.25;
} else {
s(i)= ss;
}
}
error, "no convergence after "+pr1(2^i)+" function evaluations";
}
func simpson(function, a, b, epsilon, notvector=)
/* DOCUMENT integral= simpson(function, a, b)
or integral= simpson(function, a, b, epsilon)
returns the integral of FUNCTION(x) from A to B. If EPSILON is
given, Simpson's rule is refined until that fractional accuracy
is obtained. EPSILON defaults to 1.e-6.
If the notvector= keyword is supplied and non-zero, then FUNCTION
may not be called with a list of x values to return a list of
results. By default, FUNCTION is assumed to be a vector function.
If the function is very smooth, romberg may work better.
SEE ALSO: romberg, max_doublings
*/
{
if (!epsilon || epsilon<0.) epsilon= 1.e-6;
a= double(a);
b= double(b);
ost= os= -1.e-30;
for (i=1 ; i<=max_doublings ; ++i) {
st= trapezoid(function, a, b, i, notvector);
s= (4.*st - ost)/3.;
if (abs(s-os) < epsilon*abs(os)) return s;
os= s;
ost= st;
}
error, "no convergence after "+pr1(2^i)+" function evaluations";
}
local max_doublings;
/* DOCUMENT max_doublings= 20
is the maximum number of times romberg or simpson will split the
integration interval in half. Default is 20.
*/
max_doublings= 20;
func trapezoid(function, a, b, n, notvector)
{
extern _trapezoid_i, _trapezoid_s;
if (n==1) {
_trapezoid_i= 1;
return _trapezoid_s= 0.5*(b-a)*(function(a)+function(b));
}
dx= (b-a)/_trapezoid_i;
if (!notvector) {
/* conserve memory-- dont try more than 8192 points at a time */
if (_trapezoid_i <= 8192) {
s= sum(function(span(a,b,_trapezoid_i+1)(zcen)));
} else {
x= a+(indgen(8192)-0.5)*dx;
s= sum(function(x));
dx2= 8192*dx;
for (i=8193 ; i<=_trapezoid_i ; i+=8192) s+= sum(function((x+=dx2)));
}
} else {
x= a+0.5*dx;
s= function(x);
for (i=2 ; i<=_trapezoid_i ; ++i) s+= function((x+=dx));
}
_trapezoid_i*= 2;
return _trapezoid_s= 0.5*(_trapezoid_s + dx*s);
}